source('mytheme.R')
# model version
modelversion = 'gap_log_rstan'
rstanmodelPath = 'modelrlt'
modelPath = paste0(rstanmodelPath, '/models/', modelversion)
AllExpData = read.csv(paste0("../data/AllValidData.csv")) %>%filter(valid == 1)
dur <- sort(unique(AllExpData$curDur))
AllExpData$WMSize <- factor(AllExpData$WMSize, labels = c("low", "medium", "high"))
# Replace first 3 chracters "Exp" with string "Exp. "
AllExpData$Exp <- gsub("^.{0,3}", "Exp. ", AllExpData$Exp)
AllExpData[which(AllExpData$Exp == 'Exp. 4'),"Exp"] = "Exp. 4a"
AllExpData[which(AllExpData$Exp == 'Exp. 5'),"Exp"] = "Exp. 4b"
## load model Prediction results
AllDat_predY <- read.csv(paste0(modelPath, "/rlt/AllDat_predY_",modelversion,".csv"))
AllDat_predY$WMSize <- as.factor(AllDat_predY$WMSize)
levels(AllDat_predY$WMSize) = c("low", "medium", "high")
AllDat_predY$pred_Bias = AllDat_predY$mu_r - AllDat_predY$curDur
AllDat_predY$predErr = AllDat_predY$mu_r - AllDat_predY$repDur
AllDat_predY$relatErr = AllDat_predY$predErr / AllDat_predY$repDur
# rename experiments by replacing first 3 chracters "Exp" with string "Exp. "
AllDat_predY$Exp <- gsub("^.{0,3}", "Exp. ", AllDat_predY$Exp)
AllDat_predY[which(AllDat_predY$Exp == 'Exp. 4'),"Exp"] = "Exp. 4a"
AllDat_predY[which(AllDat_predY$Exp == 'Exp. 5'),"Exp"] = "Exp. 4b"
AllDat_predY$Exp = as.factor(AllDat_predY$Exp)
AllDat_predY$gap = 1
AllDat_predY[which(AllDat_predY$Gap == 2500),"gap"] = 2
AllDat_predY$gap <- factor(AllDat_predY$gap, labels = c("short", "long"))
dat_Exp4b = dplyr::group_by(AllExpData%>%filter(Exp =='Exp. 4b'), WMSize, NSub, gap) %>%
dplyr::summarize(m_WMCrr = mean(WMCrr), n = n(), se_WMCrr = sd(WMCrr)/sqrt(n-1))%>%
dplyr::group_by(gap, WMSize)%>%
dplyr::summarize( n = n(),
mean_WMCrr = mean(m_WMCrr), se_WMCrr = sd(m_WMCrr)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize', 'NSub'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'gap'. You can override using the `.groups` argument.
dat_Exp4b$gap<- factor(dat_Exp4b$gap, labels = c("short", "long"))
plt_WMCrr_Exp4b = ggplot(data = dat_Exp4b,
aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr,group =gap, color = gap))+
geom_line(aes(linetype = gap), stat = "identity",position = position_dodge(width = 0.3))+
geom_point(aes(shape = gap), stat = "identity",position = position_dodge(width = 0.3))+
geom_errorbar(width=.3, position = position_dodge(width = 0.3)) +
coord_cartesian(ylim = c(0.5, 1)) +
scale_color_manual(values=c("#999999", "#56B4E9"))+
labs("Memory load", color = "Gap", shape = "Gap", linetype = "Gap", y = "Mean accuracy in WM task in Exp. 4b") +
theme_new
plt_WMCrr_Exp4b
model <-ezANOVA(data= AllExpData%>%filter(Exp =='Exp. 4b')%>%dplyr::group_by(NSub, WMSize, gap)%>%dplyr::summarise(m_WMCrr = mean(WMCrr)), dv = m_WMCrr, wid = NSub, within = .(WMSize, gap), type = 3, detailed = T)
## `summarise()` has grouped output by 'NSub', 'WMSize'. You can override using the
## `.groups` argument.
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "gap" will be treated as numeric.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
model
## $ANOVA
## Effect DFn DFd SSn SSd F p p<.05
## 1 WMSize 2 30 5.670075e-01 0.10014543 84.9276173 4.426759e-13 *
## 2 gap 1 15 7.853375e-05 0.01091070 0.1079680 7.470119e-01
## 3 WMSize:gap 2 30 1.273336e-03 0.03673923 0.5198813 5.998504e-01
## ges
## 1 0.793236200
## 2 0.000531086
## 3 0.008541940
model$ANOVA[4] / (model$ANOVA[4] + model$ANOVA[5])
sub_WMCrr_Exp4b = AllExpData%>%filter(Exp =='Exp. 4b')%>%dplyr::group_by(NSub, WMSize, gap)%>%dplyr::summarise(mWMCrr = mean(WMCrr))%>%pivot_wider(names_from = c("WMSize", "gap"), values_from = c(mWMCrr), names_sep="_")
## `summarise()` has grouped output by 'NSub', 'WMSize'. You can override using the
## `.groups` argument.
write.csv(sub_WMCrr_Exp4b,paste0(modelPath, '/rlt/sub_WMCrr_Exp4b.csv'))
Exp.labs.1line <- c("Exp. 1 Control", "Exp. 2 Encoding", "Exp. 3 Reproduction", "Exp. 4a Both phases", "Exp. 4b Both phases\n \ \ \ \ \ \ \ \ \ \ \ \ \ with a gap")
plt_WMCrr <- ggplot(meanForPlot, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr,
group =Exp, color = Exp))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
coord_cartesian(ylim = c(0.5, 1)) +
labs(x = "Memory load", y = "Mean accuracy in WM task") +
theme_new +scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#1a9641", "#d7191c"), labels = Exp.labs.1line)+
scale_shape_manual(values = Exp.labs.1line)
plt_WMCrr_2<-ggarrange(plt_WMCrr, plt_WMCrr_Exp4b, common.legend = FALSE, ncol=2, nrow=1, widths = c(5.5,4), labels = c("a", "b"))
plt_WMCrr_2
ggsave(paste0(getwd(), "/figures/plt_WMCrr_2.png"), plt_WMCrr_2, width = 9.5, height = 4)
### generate WM correct rates
AllExpData$WMCrr <- AllExpData$TPresent == AllExpData$WMRP
m_wmp<- dplyr::group_by(AllExpData, Exp, WMSize, NSub) %>%
dplyr::summarize(m_WMCrr = mean(WMCrr), n =n(), se_WMCrr = sd(WMCrr)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
### Average Parameters
mm_Baypar <- dplyr::group_by(Bayparlist, Exp) %>%
dplyr::summarize( m_sig_s2 = mean(sig_s2),
m_sig_pr2_log = mean(sig_pr2_log),
m_ks= mean(ks),
m_kr = mean(kr),
m_ls = mean(ls),
m_ts = mean(ts),
m_mu_pr = mean(mu_pr),
m_sig_pr2 = mean(sig_pr2),
m_mu_pr_log= mean(mu_pr_log),
m_sig_mn2 = mean(sig_mn2),
n= n(),
se_sig_s2 = sd(sig_s2)/sqrt(n-1),
se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
se_ks = sd(ks)/sqrt(n-1),
se_kr = sd(kr)/sqrt(n-1),
se_ls =sd(ls)/sqrt(n-1),
se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1))
mm_Baypar
#calculate the mean reproduction biases for the five given intervals for all subjects
mpredY_sub <- dplyr::group_by(AllDat_predY, curDur, Exp, NSub, WMSize) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
sd_repDur = sd(repDur),
m_mu_r = mean(mu_r),
m_sig_r = mean(sig_r),
m_wp = mean(wp),
se_wp = sd(wp)/sqrt(n-1),
log_lik =mean(log_lik),
cv =sd_repDur/ m_repDur,
pred_cv = mean(sig_r/mu_r),
predRP_err = mean(m_mu_r-m_repDur),
predVar_err = mean(m_sig_r-sd_repDur),
predRP_rerr = mean(abs(m_mu_r-m_repDur)/m_repDur),
predVar_rerr = mean(abs(m_sig_r-sd_repDur)/sd_repDur),
predcv_err = pred_cv-cv,
predcv_rerr = mean(abs(pred_cv-cv)/cv))
## `summarise()` has grouped output by 'curDur', 'Exp', 'NSub'. You can override
## using the `.groups` argument.
write_csv(dplyr::group_by(mpredY_sub, curDur, NSub) %>%
dplyr::summarize(m_cv = mean(cv))%>%spread(curDur, m_cv), paste0(modelPath, '/rlt/m_cv.csv'))
## `summarise()` has grouped output by 'curDur'. You can override using the
## `.groups` argument.
mpredY_sub$RP_bias = mpredY_sub$m_repDur -mpredY_sub$curDur
mpredY_sub_new <- dplyr::group_by(mpredY_sub, curDur, Exp, NSub) %>%
dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(curDur, m_RP_bias)
## `summarise()` has grouped output by 'curDur', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 1'), paste0(modelPath, '/rlt/RP_Bias_exp1.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 2'), paste0(modelPath, '/rlt/RP_Bias_exp2.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 3'), paste0(modelPath, '/rlt/RP_Bias_exp3.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 4a'), paste0(modelPath, '/rlt/RP_Bias_exp4a.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 4b'), paste0(modelPath, '/rlt/RP_Bias_exp4b.csv'))
mpredY_sub_WMsize <- dplyr::group_by(mpredY_sub, WMSize, Exp, NSub) %>%
dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(WMSize, m_RP_bias)
## `summarise()` has grouped output by 'WMSize', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp. 3'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp3.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp. 4a'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4a.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp. 4b'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4b.csv'))
#### predicted data
m_predY <- mpredY_sub%>%
dplyr::group_by(Exp, curDur, WMSize) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_sd_repDur = mean(sd_repDur),
m_m_sig_r =mean(m_sig_r),
m_m_mu_r = mean(m_mu_r),
m_m_wp = mean(m_wp),
n = n(),
m_se_wp = sd(se_wp)/sqrt(n-1),
log_lik =mean(log_lik),
mpredRP_err = mean(predRP_err),
mpredVar_err = mean(predVar_err),
mpredRP_rerr = mean(predRP_rerr),
mpredVar_rerr = mean(predVar_rerr),
cv= mean(cv),
pred_cv = mean(pred_cv),
mpredcv_err = mean(predcv_err),
mpredcv_rerr = mean(predcv_rerr))
m_predY_acc = mpredY_sub%>%
dplyr::group_by(Exp) %>%
dplyr::summarize(mpred_rerr = mean(predRP_rerr)*100,
mpredVar_rerr = mean(predVar_rerr)*100,
mpredcv_rerr = mean(predcv_rerr)*100)
m_predY_acc
#extract waic and loo-cv from parameter list
m_WAIC <- dplyr::group_by(Bayparlist, Exp) %>%
dplyr::summarize(n =n(),
m_looic = mean(looic),
m_waic = mean(waic),
se_waic = sd(waic)/sqrt(n-1),
se_looic = sd(looic)/sqrt(n-1),
m_p_loo = mean(p_loo),
m_elpd_loo = mean(elpd_loo),
m_se_looic = mean(se_looic),
m_se_p_loo = mean(se_p_loo),
m_p_waic = mean(p_waic),
m_se_waic = mean(se_waic))
m_WAIC
#load test results
AllDat_newY <- read.csv(paste0(modelPath, "/rlt/AllDat_newY_",modelversion,".csv"))
AllDat_newY$WMSize <- as.factor(AllDat_newY$WMSize)
levels(AllDat_newY$WMSize) = c("low", "medium", "high")
# Replace first 3 chracters "Exp" with string "Exp. "
AllDat_newY$Exp <- gsub("^.{0,3}", "Exp. ", AllDat_newY$Exp)
AllDat_newY[which(AllDat_newY$Exp == 'Exp. 4'),"Exp"] = "Exp. 4a"
AllDat_newY[which(AllDat_newY$Exp == 'Exp. 5'),"Exp"] = "Exp. 4b"
AllDat_newY$Exp = as.factor(AllDat_newY$Exp)
AllDat_newY$gap = 1
AllDat_newY[which(AllDat_newY$Gap == 2500, AllDat_newY$Exp == "Exp. 4b"), "gap"] = 2
AllDat_newY$gap = factor(AllDat_newY$gap, labels = c("short", "long"))
Exp.labs.2lines <- c("Exp. 1\n Control", "Exp. 2\n Encoding", "Exp. 3\n Reproduction", "Exp. 4a\n Both", "Exp. 4b\n Both, with a gap")
names(Exp.labs.2lines) <- c("Exp. 1", "Exp. 2", "Exp. 3", "Exp. 4a", "Exp. 4b")
RP_bias <- ggplot(data = m_predY%>%filter(Exp !="Exp. 4b"), aes(x = curDur, y = m_m_repDur - curDur, color=as.factor(WMSize), shape = as.factor(WMSize))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data= m_newY%>%filter(Exp !="Exp. 4b"), aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
labs(x=" ", y="Reproduction bias (s)", shape ="Memory Load", color = "Memory Load")+theme_new+
colorSet3+
theme(legend.position = "top")+ylim(-0.5, 0.4)
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias.png"), RP_bias, width = 6, height = 3)
RP_bias
RP_bias_Exp4b <- ggplot(data = AllDat_predY %>% filter(Exp == 'Exp. 4b')%>%
dplyr::group_by(Exp, curDur, WMSize, gap) %>%
dplyr::summarize(m_mu_r = mean(mu_r), m_repDur = mean(repDur), n= n(), se_mu_r = sd(mu_r)/sqrt(n-1)), aes(x = curDur, y = m_repDur - curDur, group = interaction(gap, WMSize), color=as.factor(WMSize), shape = gap)) +
geom_point(size=2, alpha = 0.5)+
geom_line(data= AllDat_newY %>% filter(Exp == 'Exp. 4b') %>%
dplyr::group_by(Exp, curDur, WMSize, gap) %>%
dplyr::summarize(m_mu_r = mean(mu_r), n= n(), se_mu_r = sd(mu_r)/sqrt(n-1)), aes(x=curDur, y=m_mu_r-curDur, group = interaction(WMSize, gap), linetype = gap, color=WMSize)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
labs(x=" ", y="Reproduction bias (s)", shape ="Gap", linetype = "Gap")+theme_new+ guides(color = "none")+
colorSet3+
theme(legend.position = "top")+ylim(-0.5, 0.4)
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
RP_bias_Exp4b
## Figures in the MS
RP_bias_all<-ggarrange(RP_bias, RP_bias_Exp4b, common.legend = FALSE, ncol=2, nrow=1, widths = c(4,1.4), labels = c("a", "b"))
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_all.png"), RP_bias_all, width = 9, height = 4.5)
RP_bias_all
RP_bias_1_4b <- ggplot(data = AllDat_predY %>%
dplyr::group_by(Exp, curDur, WMSize, gap) %>%
dplyr::summarize(m_mu_r = mean(mu_r), m_repDur = mean(repDur), n= n(), se_mu_r = sd(mu_r)/sqrt(n-1)), aes(x = curDur, y = m_repDur - curDur, group = interaction(gap, WMSize), color=as.factor(WMSize), shape = as.factor(WMSize))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data= AllDat_newY %>%
dplyr::group_by(Exp, curDur, WMSize, gap) %>%
dplyr::summarize(m_mu_r = mean(mu_r), n= n(), se_mu_r = sd(mu_r)/sqrt(n-1)), aes(x=curDur, y=m_mu_r-curDur, group = interaction(WMSize, gap), linetype = gap, color=WMSize)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
labs(x=" ", y="Reproduction bias (s)", shape ="Memory Load", linetype = "Gap", color = "Memory Load")+theme_new+
colorSet3+
theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_1_4b.png"), RP_bias_1_4b, width = 9, height = 4)
RP_bias_1_4b
AllExpData[which(AllExpData$Exp != "Exp. 4b"), "gap"] = 1
AllExpData$gap = factor(AllExpData$gap, labels = c("short", "long"))
m_obsRP = AllExpData %>% dplyr::group_by(Exp, curDur, WMSize, NSub, gap) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
se_repDur = sd(repDur)/sqrt(n-1)) %>%
dplyr::group_by(Exp, curDur, WMSize, gap) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_se_repDur = mean(se_repDur))
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize', 'NSub'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
RP_bias_obs_gap <- ggplot(data = m_obsRP, aes(x = curDur, y = m_m_repDur-curDur,shape = as.factor(WMSize), color=as.factor(WMSize), linetype = factor(gap))) +
geom_point()+
geom_line()+
geom_errorbar(width=.1, aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur, linetype = gap)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape="Memory Load", color = "Memory Load", linetype = "Gap")+
theme_new+colorSet3+
scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs_gap.png"), RP_bias_obs_gap, width = 9, height = 4)
RP_bias_obs_gap
RP_bias_obs <- ggplot(data = m_obsRP%>%filter(Exp != 'Exp. 4b'), aes(x = curDur, y = m_m_repDur-curDur,shape = as.factor(WMSize), color=as.factor(WMSize))) +
geom_point()+
geom_line()+
geom_errorbar(width=.1, aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape="Memory Load", color = "Memory Load")+
theme_new+colorSet3+
scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
RP_bias_obs
RP_bias_obs_4b <- ggplot(data = m_obsRP%>%filter(Exp == 'Exp. 4b'), aes(x = curDur, y = m_m_repDur-curDur,shape = as.factor(WMSize), color=as.factor(WMSize), linetype = factor(gap))) +
geom_point()+
geom_line()+
geom_errorbar(width=.1, aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur, linetype = gap)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
labs(x="Sample intervals (s)", y="Reproduction bias(s)", linetype = "Gap")+
theme_new+colorSet3+
scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")+ylim(-0.5, 0.2)+guides(color = "none") +guides(shape = "none")
RP_bias_obs_4b
## Figures in the MS
RP_bias_obs_all<-ggarrange(RP_bias_obs, RP_bias_obs_4b, common.legend = FALSE, ncol=2, nrow=1, widths = c(4,1.4), labels = c("a", "b"))
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs_all.png"), RP_bias_obs_all, width = 9, height = 4.5)
RP_bias_obs_all
plt_wp <- ggplot(data = AllDat_predY %>%dplyr::group_by(NSub, Exp, WMSize) %>% dplyr::summarise(m_wp = mean(wp)) %>%dplyr::group_by(Exp, WMSize) %>% dplyr::summarise(mm_wp = mean(m_wp), n= n(), m_se_wp = sd(m_wp)/sqrt(n-1)), aes(Exp, mm_wp, ymin = mm_wp - m_se_wp, ymax = mm_wp + m_se_wp, group =interaction(Exp, WMSize), color = WMSize, shape = as.factor(WMSize)))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
colorSet5+
labs(x = "", y = TeX("Weight of the prior $w_p$"), color = 'Memory Load', shape = 'Memory Load') + #scale_x_discrete(labels= Exp.labs.2lines)+
theme_new + theme(legend.position="top")
## `summarise()` has grouped output by 'NSub', 'Exp'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups` argument.
plt_wp
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_wp.png"), plt_wp, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
m_wp = AllDat_predY %>%dplyr::group_by(NSub, Exp, WMSize, gap) %>% dplyr::summarise(m_wp = mean(wp)) %>%dplyr::group_by(Exp, WMSize, gap) %>% dplyr::summarise(mm_wp = mean(m_wp), n= n(), m_se_wp = sd(m_wp)/sqrt(n-1))
## `summarise()` has grouped output by 'NSub', 'Exp', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the `.groups` argument.
plt_wp_gap <- ggplot(data = m_wp, aes(Exp, mm_wp, ymin = mm_wp - m_se_wp, ymax = mm_wp + m_se_wp, group =interaction(Exp, WMSize, gap), color = WMSize, linetype = gap, shape = WMSize))+
geom_line(stat = "identity", position = position_dodge(width = 0.3))+
geom_point(stat = "identity",position = position_dodge(width = 0.3))+
geom_errorbar(width=.3, position = position_dodge(width = .3)) +
colorSet5+
labs(x = "", y = TeX("Weight of the prior $w_p$"), color = 'Memory Load', shape = 'Memory Load', linetype= 'Gap') +
theme_new + theme(legend.position="top")
plt_wp_gap
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
sub_wp = AllDat_predY %>%dplyr::group_by(NSub, Exp, gap) %>% dplyr::summarise(wp = mean(wp))
## `summarise()` has grouped output by 'NSub', 'Exp'. You can override using the
## `.groups` argument.
Bayparlist$gap = 'short'
Bayparlist_4b = Bayparlist%>%filter(Exp == 'Exp. 4b')
Bayparlist_4b$gap = 'long'
Bayparlist_all = rbind(Bayparlist, Bayparlist_4b)
# Indifference Point based on Equation 10 in manuscript
Eq_Inp_list <- dplyr::left_join(Bayparlist_all, sub_wp) %>%
dplyr::group_by(NSub, Exp) %>%
mutate(sig_post = sig_pr2_log * wp) %>%
mutate(InP1_log = (kr*log(1)-(1-wp)*ks*log(1)+ sig_post/2 +wp*mu_pr_log)/wp)%>%
mutate(InP1 = exp(InP1_log)) %>%
mutate(InP3_log = (kr*log(3)-(1-wp)*ks*log(3)+ sig_post/2 +wp*mu_pr_log)/wp)%>%
mutate(InP3 = exp(InP3_log)) %>%
mutate(InP5_log = (kr*log(5)-(1-wp)*ks*log(5)+ sig_post/2 +wp*mu_pr_log)/wp)%>%
mutate(InP5 = exp(InP5_log)) %>%
dplyr::select(NSub, Exp,gap, InP1, InP3, InP5)
## Joining, by = c("NSub", "Exp", "gap")
write.csv(Eq_Inp_list, paste0(modelPath, '/rlt/Eq_Inp_list.csv'))
m_Eq_Inp_list = Eq_Inp_list %>%
pivot_longer(cols = starts_with("InP"), names_to = "WMSize", names_prefix = "InP", values_to = "InP")
m_Eq_Inp_list$gap = factor(m_Eq_Inp_list$gap, levels = c("short", "long"))
m_Eq_Inp_list$WMSize <- factor(m_Eq_Inp_list$WMSize, labels = c("low", "medium", "high"))
plt_Inp_eq = m_Eq_Inp_list%>%dplyr::group_by(Exp, WMSize, gap)%>%
dplyr::summarise(m_inP= mean(InP), n =n(), se_inP = sd(InP)/sqrt(n()-1)) %>%ggplot(aes(x= Exp, y=m_inP, color = WMSize, shape = WMSize, linetype = factor(gap)))+
geom_line(stat = "identity",position = position_dodge(width = 0.3))+
geom_point(stat = "identity",position = position_dodge(width = 0.3))+
geom_errorbar(width=.3, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.3)) +theme_new+
labs(colour = "Memory Load", shape = "Memory Load", linetype = "Gap")+colorSet3+
xlab(' ')+ylab("Indifference point (s)")+
theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
plt_Inp_eq
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Inp_eq.svg"), plt_Inp_eq, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
dat = Eq_Inp_list%>%filter(Exp %in%c('Exp. 2', 'Exp. 3'))
dat$high_low_2 = 0.5 * (dat$InP5 - dat$InP1)
dat1 = dat %>% dplyr::select(c("Exp","NSub","high_low_2")) %>%pivot_wider(names_from =c('Exp'), values_from = c(high_low_2))
dat1$Exp2_3 = dat1$`Exp. 2` - dat1$`Exp. 3`
mean(dat1$Exp2_3)
## [1] -0.1490995
sd(dat1$Exp2_3)/sqrt(15)
## [1] 0.03349644
mean(dat1$`Exp. 3`)
## [1] 0.07411996
sd(dat1$`Exp. 3`)/sqrt(15)
## [1] 0.02502439
AllDat_newY$predErr = AllDat_newY$mu_r - AllDat_newY$curDur
temp_newY <- AllDat_newY %>% filter(!(Exp == 'Exp. 4a' & Gap == 2500)) %>% select(Exp, WMSize, NSub, predErr, curDur, gap)
InP_curve<- temp_newY %>% dplyr::group_by(Exp, WMSize, NSub, gap)%>%
dplyr::summarise(minErr = min(abs(predErr)), idx = which.min(abs(predErr)))
## `summarise()` has grouped output by 'Exp', 'WMSize', 'NSub'. You can override
## using the `.groups` argument.
InP_curve$InP_curve = temp_newY[InP_curve$idx,]$curDur
InP_curve$y = temp_newY[InP_curve$idx,]$predErr + temp_newY[InP_curve$idx,]$curDur
#plot indifference points (the intersections of the Prediction curve with the diagonal)
plt_InP_curve<- ggplot(data = InP_curve%>%dplyr::group_by(Exp, WMSize, gap)%>% dplyr::summarise(m_InP = mean(InP_curve), se_InP = sd(InP_curve)/sqrt(n()-1)), aes(x= Exp, y=m_InP, color = as.factor(WMSize), linetype = factor(gap), shape = as.factor(WMSize)))+
geom_line(stat = "identity",position = position_dodge(width = 0.3))+
geom_point(stat = "identity",position = position_dodge(width = 0.3))+
geom_errorbar(width=.3, aes(ymin = m_InP - se_InP, ymax = m_InP + se_InP), position = position_dodge(width = 0.3)) +theme_new+
labs(colour = "Memory Load", shape ="Memory Load")+colorSet3+
xlab(' ')+ylab("indifference point (s)")+
theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_curve.png"), plt_InP_curve, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP_curve
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
#Observed Indifference Point for Exp. 4b
obs_model <- function(df) {
lm(repDur ~ curDur, data = df)
}
#Observed Indifference Point
obs_Inp_list <- AllDat_predY %>%
dplyr::group_by(NSub, Exp, WMSize, gap) %>% nest() %>%
mutate(model = map(data, obs_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
dplyr::select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, slope = curDur) # rename columns
## Warning: The `.drop` argument of `unnest()` is deprecated as of tidyr 1.0.0.
## All list-columns are now preserved.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
obs_Inp_list$model = NULL
obs_Inp_list$data = NULL
obs_Inp_list$inP = obs_Inp_list$Intercept /(1-obs_Inp_list$slope)
obs_Inp_list$slope = -1 * obs_Inp_list$slope
obs_Inp_list_no_gap <- AllDat_predY %>%
dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>%
mutate(model = map(data, obs_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
dplyr::select(-std.error,-statistic, -p.value) %>%
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, slope = curDur) # rename columns
obs_Inp_list_no_gap$model = NULL
obs_Inp_list_no_gap$data = NULL
obs_Inp_list_no_gap$inP = obs_Inp_list_no_gap$Intercept /(1-obs_Inp_list_no_gap$slope)
obs_Inp_list_no_gap$slope = -1 * obs_Inp_list_no_gap$slope
ezANOVA(data = obs_Inp_list%>%filter(Exp =='Exp. 4b'), dv= inP, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 gap 1 15 4.065468 0.062041074 0.019543220
## 3 WMSize 2 30 8.762649 0.001006806 * 0.078253848
## 4 gap:WMSize 2 30 0.584382 0.563670241 0.003061788
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.4515805 0.003829536 *
## 4 gap:WMSize 0.5210587 0.010428128 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.6458198 0.004991546 * 0.6808339 0.004255437 *
## 4 gap:WMSize 0.6761593 0.502620038 0.7194299 0.512209537
# plot the observed indifference points and slopes of RP
plt_obs_InP_slope_err<- ggplot(data = obs_Inp_list %>% group_by(Exp, WMSize)%>%
dplyr::summarise(n=n(),
m_inP = mean(inP),
se_inP = sd(inP)/sqrt(n-1),
m_slope = mean(slope),
se_slope = sd(slope)/sqrt(n-1)), aes(x= m_slope, y=m_inP, color = WMSize, shape = WMSize))+
geom_line(stat = "identity")+
geom_point(stat = "identity")+
geom_errorbar(width = 0.02, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP)) +
geom_errorbarh(height =0.02, aes(xmin = m_slope - se_slope, xmax = m_slope + se_slope)) +
theme_new+
labs(colour = "Memory Load", shape = "Memory Load")+colorSet3+
facet_grid(~Exp)+
xlab('slope of reproduction bias')+ylab("indifference point (s)")+
theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_obs_InP_slope_err.png"), plt_obs_InP_slope_err, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_obs_InP_slope_err
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
#Predicated Indifference Point for Exp.4b
pred_model <- function(df) {
lm(mu_r ~ curDur, data = df)
}
pred_Inp_list <- AllDat_predY %>%
dplyr::group_by(NSub, Exp, WMSize, gap) %>% nest() %>%
mutate(model = map(data, pred_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
dplyr::select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, pred_slope = curDur) # rename columns
pred_Inp_list$model = NULL
pred_Inp_list$data = NULL
pred_Inp_list$pred_inP = pred_Inp_list$Intercept /(1-pred_Inp_list$pred_slope)
pred_Inp_list$pred_slope = -1 * pred_Inp_list$pred_slope
pred_Inp_slope_no_gap <- AllDat_predY %>%
dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>%
mutate(model = map(data, pred_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
dplyr::select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, pred_slope = curDur) # rename columns
pred_Inp_slope_no_gap$model = NULL
pred_Inp_slope_no_gap$data = NULL
pred_Inp_slope_no_gap$pred_inP = pred_Inp_slope_no_gap$Intercept /(1-pred_Inp_slope_no_gap$pred_slope)
pred_Inp_slope_no_gap$pred_slope = -1 * pred_Inp_slope_no_gap$pred_slope
m_pred_Inp_slope_no_gap = pred_Inp_slope_no_gap %>% group_by(Exp, WMSize)%>%
dplyr::summarise(n=n(),
m_Intercept = mean(Intercept),
se_Intercept= sd(Intercept)/sqrt(n-1),
m_pred_inP = mean(pred_inP),
se_pred_inP = sd(pred_inP)/sqrt(n-1),
m_pred_slope = mean(pred_slope),
se_pred_slope = sd(pred_slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
# plot the observed indifference points and slopes of RP
plt_pred_InP_slope_err<- ggplot(data = m_pred_Inp_slope_no_gap, aes(x= m_pred_slope, y=m_pred_inP, color = WMSize, shape = WMSize))+
geom_line(stat = "identity")+geom_point(stat = "identity")+
geom_errorbar(width = 0.02, aes(ymin = m_pred_inP - se_pred_inP, ymax = m_pred_inP + se_pred_inP)) +
geom_errorbarh(height =0.02, aes(xmin = m_pred_slope - se_pred_slope, xmax = m_pred_slope + se_pred_slope)) +
geom_point(data = obs_Inp_list%>% group_by(Exp, WMSize)%>%
dplyr::summarise(n=n(),
m_inP = mean(inP),
se_inP = sd(inP)/sqrt(n-1),
m_slope = mean(slope),
se_slope = sd(slope)/sqrt(n-1)), aes(x= m_slope, y =m_inP, color = WMSize, shape = WMSize))+
theme_new+
labs(colour = "Memory Load", shape = "Memory Load")+colorSet3+
facet_grid(~Exp)+
xlab('slope of reproduction')+ylab("indifference point (s)")+
theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_pred_InP_slope_err.png"), plt_pred_InP_slope_err, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_pred_InP_slope_err
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
InP_obs<- ggplot(data = obs_Inp_list_no_gap %>%dplyr::group_by(WMSize, Exp) %>%dplyr::summarise(m_inP = mean(inP), se_inP = sd(inP)/sqrt(n()-1)), aes(x= Exp, y=m_inP, color = WMSize, shape = WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.3))+
geom_point(stat = "identity",position = position_dodge(width = 0.3))+
geom_errorbar(width=.3, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.3)) +theme_new+
labs(colour = "Memory Load", shape = "Memory Load")+colorSet3+
xlab(' ')+ylab("observed indifference point (s)")+
theme(legend.position = "top")
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/InP_obs.png"), InP_obs, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
InP_obs
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
InP_pred<- ggplot(data = m_pred_Inp_slope_no_gap, aes(x= Exp, y=m_pred_inP, color = WMSize, shape = WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.3))+
geom_point(stat = "identity",position = position_dodge(width = 0.3))+
geom_errorbar(width=.3, aes(ymin = m_pred_inP - se_pred_inP, ymax = m_pred_inP + se_pred_inP), position = position_dodge(width = 0.3)) +theme_new+
labs(colour = "Memory Load", shape = "Memory Load")+colorSet3+ #scale_x_discrete(labels= Exp.labs.2lines)+
xlab(' ')+ylab("indifference point (s)")+
theme(legend.position = "top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/InP_pred.png"), InP_pred, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
InP_pred
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
### Calculate predication error
Inp_list_no_gap = left_join(obs_Inp_list_no_gap, pred_Inp_slope_no_gap, by = c("NSub", "Exp", "WMSize"))
Inp_list_no_gap$InP_err = Inp_list_no_gap$pred_inP -Inp_list_no_gap$inP
Inp_list_no_gap$InP_rerr = 100*Inp_list_no_gap$InP_err/ Inp_list_no_gap$inP
Inp_list_no_gap$slope_err = Inp_list_no_gap$pred_slope - Inp_list_no_gap$slope
Inp_list_no_gap$slope_rerr = 100* Inp_list_no_gap$slope_err/Inp_list_no_gap$slope
m_Inp_list_no_gap = Inp_list_no_gap %>% dplyr::group_by(Exp) %>% dplyr::summarise(m_InP_rerr = mean(InP_rerr), m_slope_rerr = mean(slope_rerr), m_InP_rerr_abs = mean(abs(InP_rerr)), m_slope_rerr_abs = mean(abs(slope_rerr)))
m_Inp_list_no_gap$InP_auc = 100- m_Inp_list_no_gap$m_InP_rerr_abs
m_Inp_list_no_gap$slope_auc = 100- m_Inp_list_no_gap$m_slope_rerr_abs
x
#plot the predicated indifference points and slope of predicated RP
plt_pred_InP_slope_err<- ggplot(data = obs_Inp_list_no_gap%>% dplyr::group_by(Exp, WMSize)%>%
dplyr::summarise(n=n(),
m_inP = mean(inP),
se_inP = sd(inP)/sqrt(n-1),
m_slope = mean(slope),
se_slope = sd(slope)/sqrt(n-1)), aes(x= m_slope, y=m_inP, shape = WMSize, color = WMSize))+
geom_line(stat = "identity")+
geom_point(stat = "identity")+
geom_errorbar(width = 0.02, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP)) +
geom_errorbarh(height =0.02, aes(xmin = m_slope - se_slope, xmax = m_slope + se_slope)) +
theme_new+
labs(colour = "Memory Load", shape = "Memory Load") +colorSet3+
facet_grid(~Exp)+
xlab('slope of reproduction')+ylab("indifference point (s)")+
theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_pred_InP_slope_err.png"), plt_pred_InP_slope_err, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_pred_InP_slope_err
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## Figures in the MS
fig3<-ggarrange(RP_bias, plt_pred_InP_slope_err, common.legend = TRUE, ncol=1, nrow=2, labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig3.png"), fig3, width = 6, height = 5)
fig3
## combine InP and wp
fig_wp_InP_gap<-ggarrange(plt_wp_gap, plt_Inp_eq, common.legend = TRUE, ncol=2, nrow=1, labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_wp_InP_gap.png"), fig_wp_InP_gap, width = 6, height = 3)
fig_wp_InP_gap
## combine InP and wp
fig4<-ggarrange(plt_wp, InP_pred, common.legend = TRUE, ncol=2, nrow=1, labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig4.png"), fig4, width = 6, height = 3)
fig4
m_predErr_sub<- mpredY_sub%>%
dplyr::group_by(Exp, WMSize, NSub) %>% dplyr::summarise(
mpredRP_err=mean(predRP_err),
mpredVar_err=mean(predVar_err),
mpredcv_err = mean(predcv_err),
mpredRP_rerr = mean(predRP_rerr),
mpredVar_rerr = mean(predVar_rerr),
mpredcv_rerr = mean(predcv_rerr))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
m_predErr<- m_predY%>%
dplyr::group_by(Exp, WMSize) %>% dplyr::summarise(
mmpredcv_err = mean(mpredcv_err),
mmpredRP_err=mean(mpredRP_err),
mmpredVar_err=mean(mpredVar_err),
mmpredRP_rerr = mean(mpredRP_rerr),
mmpredVar_rerr = mean(mpredVar_rerr),
mmpredcv_rerr = mean(mpredcv_rerr))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predErr
m_predErr_sub$model = 'logarithmic'
m_predErr$model = 'logarithmic'
linear_model = 'gap_linear_rstan'
m_predErr_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_", linear_model, ".csv"))
m_predErr_linear$X = NULL
m_predErr_sub_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_sub_", linear_model, ".csv"))
m_predErr_sub_linear$X = NULL
m_predErr_sub_all = rbind(m_predErr_sub, m_predErr_sub_linear)
m_predErr_all = rbind(m_predErr, m_predErr_linear)
m_predErr_all$WMSize = as.factor(m_predErr_all$WMSize)
levels(m_predErr_all$WMSize) = c("low", "medium", "high")
temp = m_predErr_all %>% filter(model == 'logarithmic') %>%summarise(abs_mmpredcv_err = abs(mmpredcv_err))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_Err_CV_all = ggplot(m_predErr_all, aes(abs(mmpredRP_err), abs(mmpredcv_err), group = interaction(model, Exp, WMSize), color = model, shape = Exp, size = WMSize)) +
geom_point() +
geom_hline(yintercept = round(max(temp$abs_mmpredcv_err), 4)+0.0005, linetype='dashed')+
xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in CV')+colorSet3+
scale_shape_manual(values = c(6, 7, 16, 17,8)) +
theme_new+
theme(legend.position = 'top')+
labs(size = 'Memory Load')+
guides(colour = guide_legend(order = 1, nrow=2,byrow=TRUE),
shape = guide_legend(order =2, nrow=2,byrow=TRUE),
size = guide_legend(order = 3, nrow=3,byrow=TRUE))
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Err_CV_all.png"), plt_Err_CV_all, width = 7, height = 4)
## Warning: Using size for a discrete variable is not advised.
plt_Err_CV_all
## Warning: Using size for a discrete variable is not advised.
m_predY_acc = m_predErr_sub_all%>%
dplyr::group_by(Exp, model) %>%
dplyr::summarize(mmpredRP_rerr = mean(mpredRP_rerr)*100,
mmpredVar_rerr = mean(mpredVar_rerr)*100,
mmpredcv_rerr = mean(mpredcv_rerr)*100,
mmpredRP_acc = (1-mean(mpredRP_rerr))*100,
mmpredVar_acc = (1-mean(mpredVar_rerr))*100,
mmpredCV_acc = (1-mean(mpredcv_rerr))*100)
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predY_acc